Curso I – SISSA

Archivos raster y datos espaciales

Objetivo y contenido

Objetivo y contenido

Objetivo: aprender a manipular conjuntos de datos grillados y datos vectoriales mediante la aplicación de funciones en R.

  • Importar y trazar datos vectoriales
  • Extraer información espacial específica utilizando datos vectoriales
  • Importar conjuntos de datos grillados
  • Manipular conjuntos de datos grillados
  • Exportar conjuntos de datos grillados en formatos diferentes
  • Manipular valores en datos grillados

Importar y visualizar datos vectoriales

Datos vectoriales

Los datos vectoriales almacenan la ubicación, forma y atributos de características geográficas en una de las siguientes formas:

  • Puntos
  • Líneas
  • Polígonos

Existen diferentes paquetes que se pueden utilizar para cargar y manipular datos vectoriales (shapefiles) en R. En esta presentación, utilizaremos los paquetes sf y terra. El paquete tibble también puede ser instalado ya que es utilizado por sf.

Datos vectoriales

Para instalar estos paquetes podemos utilizar la función install.packages.

install.packages("sf")
install.packages("tibble")
install.packages("terra")

Una vez instalados, podemos cargarlos con la función library.

library(sf)
library(tibble)
library(terra)

Datos vectoriales

Para leer datos vectoriales en R, podemos utilizar ambos paquetes. Para utilizar el paquete sf podemos aplicar la función read_sf.

countries <- read_sf("C:/User/Data/L3_Rasters_and_Shapefiles_Data/Countries")

Si la carpeta contiene más de un shapefile, hay que especificar la ruta del archivo:

countries <- read_sf("C:/User/Data/L3_Rasters_and_Shapefiles_Data/Countries/Nile_Basin_Countries_GAUL_2014_2.shp")

Datos vectoriales

Si escribimos el nombre del objeto en la consola:

countries
## Simple feature collection with 15 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 12.19545 ymin: -13.45568 xmax: 47.98618 ymax: 31.66734
## Geodetic CRS:  WGS 84
## # A tibble: 15 × 10
##    STATUS           DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR Shape_Leng
##    <chr>            <chr>         <int> <chr>         <int>     <int>      <dbl>
##  1 Member State     NO              205 Rwanda         1000      3000       8.13
##  2 Member State     NO              133 Kenya          1000      3000      48.5 
##  3 Member State     NO               74 South Su…      2011      3000      46.9 
##  4 Member State     NO              257 United R…      1000      3000      83.0 
##  5 Member State     NO              253 Uganda         1000      3000      23.2 
##  6 Sovereignty uns… YES           61013 Ilemi tr…      1000      3000       2.84
##  7 Member State     NO               79 Ethiopia       1000      3000      50.4 
##  8 Member State     NO               77 Eritrea        1000      3000      50.0 
##  9 Member State     NO               43 Burundi        1000      3000       9.12
## 10 Member State     NO               68 Democrat…      1000      3000      99.0 
## 11 Sovereignty uns… YES           40760 Hala'ib …      1000      3000       8.51
## 12 Sovereignty uns… YES           40762 Ma'tan a…      1000      3000       2.06
## 13 Member State     NO                6 Sudan          2011      3000      81.9 
## 14 Member State     NO            40765 Egypt          1000      3000      61.3 
## 15 Sovereignty uns… YES             102 Abyei          2011      3000       3.61
## # ℹ 3 more variables: Shape_Area <dbl>, Name_label <chr>,
## #   geometry <MULTIPOLYGON [°]>

Datos vectoriales

Vamos a graficar el objeto countries:

plot(countries)

Datos vectoriales

La función names devolverá los nombres de todos los atributos almacenados.

names(countries)
##  [1] "STATUS"     "DISP_AREA"  "ADM0_CODE"  "ADM0_NAME"  "STR0_YEAR" 
##  [6] "EXP0_YEAR"  "Shape_Leng" "Shape_Area" "Name_label" "geometry"

Simple Features en R

Los objetos sf son objetos que cuentan con una columna especial de “geometría”:

st_geometry(countries)
## Geometry set for 15 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 12.19545 ymin: -13.45568 xmax: 47.98618 ymax: 31.66734
## Geodetic CRS:  WGS 84
## First 5 geometries:
## MULTIPOLYGON (((30.47352 -1.057411, 30.47035 -1...
## MULTIPOLYGON (((39.37656 -4.721247, 39.37412 -4...
## MULTIPOLYGON (((33.02315 12.22673, 33.25586 12....
## MULTIPOLYGON (((40.21717 -10.29314, 40.21717 -1...
## MULTIPOLYGON (((34.09062 3.878785, 34.09114 3.8...

Simple Features en R

Ya que las geometrías no contienen valorer únicos, se almacenan en una columna-lista, con cada elemento en la lista conteniendo la geometría de ese elemento. Las tres clases usadas para representar Simple Features son:

  • sf, la tabla (data.frame) con los atributos y geometrías
  • sfc, la columna-lista con las geometrías para cada elemento
  • sfg, la geometría de cada elemento individual

Simple Features en R

Si queremos graficar solo la geometría:

countries_geom <- st_geometry(countries)
plot(countries_geom, axes = TRUE)

Simple Features en R

También podemos definir los límites para el eje x y el eje y del gráfico.

countries_geom <- st_geometry(countries)
plot(st_geometry(countries), axes = TRUE,
xlim = c(20, 40), ylim = c(1, 20))

Simple Features en R

Si queremos graficar un solo atributo:

plot(countries["Shape_Area"], axes=TRUE)

Existen muchas más posibilidades: https://r-spatial.github.io/sf/articles/sf5.html

Extracción de información espacial específica usando datos vectoriales

Extrayendo Polígonos (o Features) de Shapefiles

Cada fila contiene una única característica:

Extrayendo Polígonos (o Features) de Shapefiles

Para obtener una simple feature que contenga solo Etiopía:

ethiopia <- countries[7, ]
print(ethiopia)
## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 32.99994 ymin: 3.401109 xmax: 47.98618 ymax: 14.89416
## Geodetic CRS:  WGS 84
## # A tibble: 1 × 10
##   STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR Shape_Leng Shape_Area
##   <chr>  <chr>         <int> <chr>         <int>     <int>      <dbl>      <dbl>
## 1 Membe… NO               79 Ethiopia       1000      3000       50.4       92.9
## # ℹ 2 more variables: Name_label <chr>, geometry <MULTIPOLYGON [°]>

Extrayendo Polígonos (o Features) de Shapefiles

plot(ethiopia)

Exportando datos vectoriales

Para guardar la(s) característica(s) extraída(s) en una carpeta específica del sistema, podemos usar la función write_sf.

write_sf(ethiopia, "C:/User/Ethopida_outline/Ethiopia_outline.shp")

Puede guardar en otros formatos con st_write. Consulte el argumento driver en la ayuda.

Manipulando atributos

Podemos deshacernos de la geometría con la función st_drop_geometry.

Calculando nuevos atributos (por ejemplo, mm/m2 al total por país):

countries <- read_sf("C:/User/Data/L3_Rasters_and_Shapefiles_Data/Countries")
countries["precip_mm"]    <- runif(nrow(countries))
countries_df              <- st_drop_geometry(countries)
countries["precip_total"] <- countries_df["precip_mm"] * countries_df["Shape_Area"]

Manipulando atributos

Alternativamente, podemos usar dplyr:

countries %>%
mutate(precip_ml = runif(nrow(countries)),
precip_total = precip_ml * Shape_Area)
## Simple feature collection with 15 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 12.19545 ymin: -13.45568 xmax: 47.98618 ymax: 31.66734
## Geodetic CRS:  WGS 84
## # A tibble: 15 × 13
##    STATUS           DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR Shape_Leng
##  * <chr>            <chr>         <int> <chr>         <int>     <int>      <dbl>
##  1 Member State     NO              205 Rwanda         1000      3000       8.13
##  2 Member State     NO              133 Kenya          1000      3000      48.5 
##  3 Member State     NO               74 South Su…      2011      3000      46.9 
##  4 Member State     NO              257 United R…      1000      3000      83.0 
##  5 Member State     NO              253 Uganda         1000      3000      23.2 
##  6 Sovereignty uns… YES           61013 Ilemi tr…      1000      3000       2.84
##  7 Member State     NO               79 Ethiopia       1000      3000      50.4 
##  8 Member State     NO               77 Eritrea        1000      3000      50.0 
##  9 Member State     NO               43 Burundi        1000      3000       9.12
## 10 Member State     NO               68 Democrat…      1000      3000      99.0 
## 11 Sovereignty uns… YES           40760 Hala'ib …      1000      3000       8.51
## 12 Sovereignty uns… YES           40762 Ma'tan a…      1000      3000       2.06
## 13 Member State     NO                6 Sudan          2011      3000      81.9 
## 14 Member State     NO            40765 Egypt          1000      3000      61.3 
## 15 Sovereignty uns… YES             102 Abyei          2011      3000       3.61
## # ℹ 6 more variables: Shape_Area <dbl>, Name_label <chr>,
## #   geometry <MULTIPOLYGON [°]>, precip_mm <dbl>, precip_total <dbl>,
## #   precip_ml <dbl>

Lectura de Datos vectoriales con el paquete terra

Cuando queremos manipular rásters y datos vectoriales simultaneamente, el paquete terra es muy útil. Para leer datos vectoriales con este paquete, usaremos la función vect.

countries <- vect("C:/User/Data/L3_Rasters_and_Shapefiles_Data/Countries/Nile_Basin_Countries_GAUL2014_2.shp")

Importando conjuntos de datos grillados

Archivos ráster

Un ráster consiste en una matriz de celdas (o píxeles), donde cada celda contiene un valor o conjunto de valores que representan cierta información.

Archivos ráster

Los ráster con datos discretos a menudo utilizan un código para representar diferentes tipos de datos.

  • Por ejemplo, estos son los valores utilizados para el producto de Cobertura de Suelo MODIS para representar la clasificación IGBP (Loveland & Belward, 1997).

Archivos ráster

Para abrir un archivo ráster GeoTIFF (o similar), se utiliza el comando rast del paquete terra.

gleam.path <- "C:/User/L3_Rasters_and_Shapefiles_Data/GLEAM_annual/GLEAM_2009.tif"
gleam      <- rast(gleam.path)

Para poder realizar operaciones ráster, todos los archivos ráster involucrados deberán tener la misma geometría:

  • Misma extensión espacial
  • Misma resolución espacial
  • Mismo origen

Archivos ráster - Sistema de coordenadas de referencia

Los metadatos del archivo ráster se pueden visualizar escribiendo el nombre del objeto donde se almacenó (en este caso, gleam).

print(gleam)
## class       : SpatRaster 
## dimensions  : 720, 1440, 1  (nrow, ncol, nlyr)
## resolution  : 0.25, 0.25  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : GLEAM_2009.tif 
## name        : GLEAM_2009 
## min value   :  -23.93291 
## max value   : 2148.19751

Archivos ráster - Sistema de coordenadas de referencia

En R, la información sobre el ráster (sistema de coordenadas, extensión espacial, número de celdas, etc.) se puede obtener o introducir fácilmente usando la función crs del paquete raster:

crs(gleam)
## [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"
crs(gleam) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"

Representación gráfica de un ráster

Para graficar un ráster usamos la función plot.

plot(gleam)

Archivos ráster - Representación gráfica de un ráster

Alternativamente, podemos usar la función click para obtener el valor de un cierto píxel.

click(gleam)

Para salir de la función click, presione Esc.

Archivos ráster

El paquete terra es capaz de trabajar eficientemente con un gran número de archivos ráster; por lo tanto, crear objetos de varias capas es extremadamente fácil.

El uso de objetos grillados de varias capas (stacks) es muy útil (por ejemplo, cuando se trabaja con rasters que involucran datos ambientales como precipitación, evapotranspiración y temperatura).

Creación de un objeto grillado de varias capas

Queremos usar estos archivos anuales de ETa para:

  • Crear una pila con todos los archivos anuales
  • Calcular un raster promedio multianual de Ea

Manipulación de conjuntos grillados

Manipulando un stack - Ejemplo

Paso 1: Listar todos los archivos dentro de la carpeta.

ssebop.path <- "C:/User/Example/SSEBop_annual"
ssebop.files <- list.files(ssebop.path, full.names = TRUE, pattern = ".tif")
head(ssebop.files)
## [1] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2003.tif"
## [2] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2004.tif"
## [3] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2005.tif"
## [4] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2006.tif"
## [5] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2007.tif"
## [6] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2008.tif"

Manipulando un stack - Ejemplo

Paso 2: Crear un stack a partir de todas las capas ráster individuales.

ssebop <- rast(ssebop.files)
plot(ssebop)

Manipulando un stack - Ejemplo

Manipulando un stack - Ejemplo

Nota: también podemos usar la función names para ver los nombres de las capas.
names(ssebop)
##  [1] "SSEBop_2003" "SSEBop_2004" "SSEBop_2005" "SSEBop_2006" "SSEBop_2007"
##  [6] "SSEBop_2008" "SSEBop_2009" "SSEBop_2010" "SSEBop_2011" "SSEBop_2012"
## [11] "SSEBop_2013" "SSEBop_2014" "SSEBop_2015" "SSEBop_2016" "SSEBop_2017"
## [16] "SSEBop_2018"

Manipulando un stack - Ejemplo

Paso 3: Calcular la media multianual.

ssebop.multiannual <- mean(ssebop)
plot(ssebop.multiannual)

Acceso a Rasters desde un Raster Stack

Actualmente tenemos el ráster stack guardado como ssebop.

print(ssebop)
## class       : SpatRaster 
## dimensions  : 8082, 7772, 16  (nrow, ncol, nlyr)
## resolution  : 0.009652, 0.009652  (x, y)
## extent      : -20.00845, 55.00689, -38.00535, 40.00211  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : SSEBop_2003.tif  
##               SSEBop_2004.tif  
##               SSEBop_2005.tif  
##               ... and 13 more source(s)
## names       : SSEBop_2003, SSEBop_2004, SSEBop_2005, SSEBop_2006, SSEBop_2007, SSEBop_2008, ... 
## min values  :           0,           0,           0,           0,           0,           0, ... 
## max values  :        3272,        3353,        3299,        3296,        3462,        3333, ...

Acceso a Rásters desde un Raster Stack

Al igual que con una lista, usamos dobles corchetes para acceder a un ráster específico (o varios rásters) desde un stack.

print(ssebop[[2]])
## class       : SpatRaster 
## dimensions  : 8082, 7772, 1  (nrow, ncol, nlyr)
## resolution  : 0.009652, 0.009652  (x, y)
## extent      : -20.00845, 55.00689, -38.00535, 40.00211  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : SSEBop_2004.tif 
## name        : SSEBop_2004 
## min value   :           0 
## max value   :        3353

Acceso a Rásters desde un Raster Stack

Al igual que con una lista, usamos dobles corchetes para acceder a un ráster específico (o varios rásters) desde un stack.

print(ssebop[[1:3]])
## class       : SpatRaster 
## dimensions  : 8082, 7772, 3  (nrow, ncol, nlyr)
## resolution  : 0.009652, 0.009652  (x, y)
## extent      : -20.00845, 55.00689, -38.00535, 40.00211  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : SSEBop_2003.tif  
##               SSEBop_2004.tif  
##               SSEBop_2005.tif  
## names       : SSEBop_2003, SSEBop_2004, SSEBop_2005 
## min values  :           0,           0,           0 
## max values  :        3272,        3353,        3299

Recortar Rásters y Extraer por Máscara

  • Recortar un ráster (función crop) extraerá la parte de un ráster dentro de un área rectangular.

  • Enmascarar un ráster (función mask), extraerá los píxeles dentro del shapefile.

Función Crop

La función crop extrae una porción de un conjunto de datos grillados según una entrada seleccionada (puede ser un archivo ráster o un archivo vectorial; es decir, SpatRaster o SpatVector, respectivamente).

Para recortar un ráster:

rwanda     <- countries[1, ]
gleam.crop <- crop(gleam, rwanda)

Función Crop

plot(gleam.crop)

Función Crop

El parámetro snap se utiliza para alinear el archivo recortado con la extensión del elemento seleccionado que se utiliza para realizar el recorte.

  • “near” (predeterminado)
  • “in”
  • “out”
gleam.crop <- crop(gleam, rwanda, snap = "out")

Función Crop

plot(gleam.crop)

Función Mask

Para enmascarar un ráster:

gleam.mask <- mask(gleam.crop, rwanda)
plot(gleam.mask)

¿Qué valores hay en mi Ráster?

Para listar todos los valores únicos en un objeto ráster, se utiliza la función values.

-Como resultado se ontiene un vector

gleam_vals <- values(gleam.mask)
gleam_vals <- as.numeric(gleam_vals)
print(gleam_vals)
##  [1]        NA        NA        NA        NA        NA  823.3318  742.6441
##  [8]  824.1913        NA        NA        NA 1063.5068 1031.8060  977.4069
## [15]  837.3028  793.3273  835.6875        NA        NA 1141.6630 1031.7339
## [22]  938.8115  869.1260  827.8035  774.6710  820.8055  905.7816        NA
## [29] 1192.4106 1073.8735  914.8868  840.9797  838.4015  800.3232  828.6646
## [36]  908.3577        NA 1179.2241 1033.6238  878.3455  838.2838  795.3875
## [43]  769.3068  844.0502  860.8229 1094.2697 1071.3359  979.1405  876.9600
## [50]  837.6367  830.4340  897.6795  830.6163  840.5397  966.6479  974.6377
## [57] 1011.7329  880.3939  856.8268        NA        NA        NA        NA
## [64]        NA        NA  957.3956  891.5947  873.2820        NA        NA
## [71]        NA        NA

Remuestrear Rásters

Para realizar operaciones ráster, la geometría entre las capas debe ser la misma.

  • Se pueden utilizar diversos métodos.

  • Aquí nos enfocaremos en dos: interpolación bilineal e interpolación por el vecino más cercano (nearest neighbour).

Remuestrear Rásters

La diapositiva anterior mostró el remuestreo a un origen diferente utilizando el mismo tamaño de celda de la grilla.

El remuestreo también se puede usar para cambiar la resolución espacial de la(s) capa(s).

Remuestrear Rásters

Cómo remuestrear rásters en R. - Método nearest neighbour (rásters con variables categóricas o continuas)

land.cover <- "C:/Users/Example/MCD12Q1_LC1_2009_001.tif"
land.cover <- rast(land.cover)
print(land.cover)
## class       : SpatRaster 
## dimensions  : 6002, 4802, 1  (nrow, ncol, nlyr)
## resolution  : 0.008329863, 0.008330556  (x, y)
## extent      : 10, 50, -15, 35  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : MCD12Q1_LC1_2009_001.tif 
## name        : MCD12Q1_LC1_2009_001 
## min value   :                    1 
## max value   :                   17

Remuestrear Rásters

Para remuestrear el ráster a la resolución espacial del objeto gleam, usaremos la función resample.

landcover.res <- resample(land.cover, gleam, method = "near")

Remuestrear Rásters

print(landcover.res)
## class       : SpatRaster 
## dimensions  : 720, 1440, 1  (nrow, ncol, nlyr)
## resolution  : 0.25, 0.25  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : MCD12Q1_LC1_2009_001 
## min value   :                    2 
## max value   :                   17

Remuestrear Rásters

plot(landcover.res, xlim = c(-25, 80), ylim = c(-15, 35))

Remuestrear Rásters

Cómo remuestrear rásters en R. Método bilineal (rásters con variables continuas).

chirps <- "C:/Users/Example/chirps_monthly2009-01"
chirps <- rast(chirps)
print(chirps)
## class       : SpatRaster 
## dimensions  : 2000, 7200, 1  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : chirps_monthly2009-01.tif 
## name        : chirps_monthly2009-01 
## min value   :                 0.000 
## max value   :              1301.122
chirps.res <- resample(chirps, gleam, method = "bilinear")

Remuestrear Rásters

print(chirps.res)
## class       : SpatRaster 
## dimensions  : 720, 1440, 1  (nrow, ncol, nlyr)
## resolution  : 0.25, 0.25  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : chirps_monthly2009-01 
## min value   :                 0.000 
## max value   :              1221.876

Remuestrear Rásters

plot(chirps.res, xlim = c(-25, 80), ylim = c(-15, 35))

Exportar conjuntos de datos grillados

Exportar conjuntos de datos grillados

Después de haber modificado un rásters, es posible que deseemos exportarlo. Esto se puede hacer con la función writeRaster (parte del paquete terra).

writeRaster (x, filename, overwrite, ...)

x: objeto ráster. filename: nombre de archivo de salida (con la ruta de archivo completa). overwrite: argumento lógico.

Exportar conjuntos de datos grillados

Utilicemos la función writeRaster para guardar el ráster CHIRPS remuestreado del último ejercicio. Guardaremos el ráster como un archivo GeoTiff.

writeRaster(chirps.res, "agregar-ruta-de-archivo" overwrite = TRUE)

Pregunta: ¿Cómo podrías hacer que el nombre de la ruta de archivo cambie según el valor de una variable?

Manipulación de datos grillados

Operaciones ráster

Una vez que los rásters tienen la misma geometría, es fácil realizar operaciones matemáticas.

  • Ejemplo: si tenemos los siguientes tres rásters:

Operaciones de ráster - valores NA values

Ahora, si tenemos estos rásters en un stack:

¿Cuál sería el resultado de las siguientes expresiones para la celda de cuadrícula resaltada:

  • mean(stacked, na.rm = TRUE)
  • mean(stacked, na.rm = FALSE)

Aislar celdas con atributos particulares

Podemos identificar píxeles con ciertos atributos. Por ejemplo, las regiones con más de 200 mm de precipitación durante enero de 2009.

plot(chirps, main = "Precipitación enero - 2009")

Aislar celdas con atributos particulares

Aislar celdas con atributos particulares

Para este propósito, estableceremos los valores mayores a 200 mm igual a 1, mientras que los valores menores o iguales a 200 mm se reemplazarán por un valor de 0.

chirps.lim <- chirps
chirps.lim[chirps.lim <= 200] <- 0
chirps.lim[chirps.lim > 200]  <- 1

Aislar celdas con atributos particulares

plot(chirps.lim, main = "Regiones con más de 200 mm para enero de 2009")

Aislar celdas con atributos particulares

plot(chirps.lim)

Aislar celdas con atributos particulares

Otro uso común de esta función es crear un ráster donde los valores son 1 o NA (de aquí en adelante, llamaremos a esto un ráster “unitario”). En este ejemplo, queremos calcular la precipitación promedio de las zónas áridas para enero de 2009 en los países del Nilo usando las zonas climáticas principales de Köppen-Geiger.

climate.zones <- "C:/Users/Example/NB_Major_CZ.tif"
climate.zones <- rast(climate.zones)

Aislar celdas con atributos particulares

plot(climate.zones,
main = "Zonas climáticas principales")

Aislar celdas con atributos particulares

Ahora Encontraremos valores únicos usando la función values en combinación con la función unique.

uniq.values <- as.numeric(unique(values(climate.zones)))
print(uniq.values)
## [1] NaN   2   0   1   3   4
Valor Clima
1 Tropical
2 Árido
3 Templado
4 Polar

Aislar celdas con atributos particulares

A continuación, establecemos los valores de la zona climática árida (valor = 2) a 1 y todo lo demás a NA.

arid <- climate.zones
arid[arid != 2] <- NA
arid[arid == 2] <- 1

Aislar celdas con atributos particulares

plot(arid, col = "red")

Aislar celdas con atributos particulares

print(arid)
## class       : SpatRaster 
## dimensions  : 10548, 8367, 1  (nrow, ncol, nlyr)
## resolution  : 0.004278282, 0.004278282  (x, y)
## extent      : 12.1931, 47.98949, -13.45713, 31.67019  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : spat_vR6JxcbrnNc7Nbj_40545.tif 
## name        : NB_Major_CZ 
## min value   :           1 
## max value   :           1
print(chirps)
## class       : SpatRaster 
## dimensions  : 2000, 7200, 1  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : chirps_monthly2009-01.tif 
## name        : chirps_monthly2009-01 
## min value   :                 0.000 
## max value   :              1301.122

Aislar celdas con atributos particulares

Vamos a redimensionar el ráster unitario a la misma resolución espacial que el producto de precipitación.

arid <- resample(arid, chirps, method = "near")
print(arid)
## class       : SpatRaster 
## dimensions  : 2000, 7200, 1  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : NB_Major_CZ 
## min value   :           1 
## max value   :           1

Aislar celdas con atributos particulares

Ahora, multiplicamos el ráster de precipitación por el ráster unitario de las zonas áridas y lo recortamos usando la capa climate.zones.

arid.precipitation <- chirps * arid
arid.precipitation <- crop(arid.precipitation, climate.zones, snap = "out")

Aislar celdas con atributos particulares

plot(arid.precipitation, main = "P para la zona árida - enero, 2009")

Aislar celdas con atributos particulares

Finalmente, queremos obtener la precipitación media para la zona climática árida.

  • Para este propósito, utilizaremos la función values para extraer los valores numéricos de la capa ráster.
  • Luego utilizaremos la función mean con el parámetro na.rm establecido en TRUE.
mean.arid <- mean(values(arid.precipitation), na.rm = TRUE)
print(mean.arid)
## [1] 6.516357

Extrayendo valores de un ráster

A veces, queremos extraer valores específicos (según las coordenadas) de un ráster. Para este propósito, utilizaremos la función extract del paquete raster. El primer paso es crear dos vectores con valores de latitud y longitud.

longitude <- c(-100, -50, 20, 100, 117)
latitude  <- c(22, -20, -1, 0, -32)
lonlat    <- cbind(longitude, latitude)
print(lonlat)
##      longitude latitude
## [1,]      -100       22
## [2,]       -50      -20
## [3,]        20       -1
## [4,]       100        0
## [5,]       117      -32

Extrayendo valores de un ráster

Ahora tenemos que convertir las coordenadas a un SpatVector.

locations <- vect(lonlat)
print(locations)
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 5, 0  (geometries, attributes)
##  extent      : -100, 117, -32, 22  (xmin, xmax, ymin, ymax)
##  coord. ref. :

Extrayendo valores de un ráster

Podemos graficar CHIRPS y las estaciones de la siguiente manera:

plot(chirps)
points(locations)

Extrayendo valores de un ráster

Para extraer los valores de las estaciones seleccionadas utilizaremos la función extract.

extraction <- terra::extract(chirps, locations)
extraction <- extraction[2:length(extraction)]
print(extraction)
##   chirps_monthly2009-01
## 1              6.863122
## 2            259.973359
## 3            183.161122
## 4            329.991304
## 5             16.177693

Extrayendo valores de un ráster

Como la función extract puede estar incluida en otro paquete, los dobles dos puntos (::) relacionan el paquete que se debe usar con la función respectiva. También podemos usar extract para extraer series temporales de un objeto stack.

¡Gracias por su atención!